<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8" />
<meta name="Generator" content="Kate, the KDE Advanced Text Editor" />
<title>tutorial9.py</title>
</head>
<body>
<pre style='color:#141312;background-color:#ffffff;'>
<i><span style='color:#008000;'>#!/usr/bin/env python</span></i>

<i><span style='color:#008000;'>&quot;&quot;&quot;********************************************************************************</span></i>
<i><span style='color:#008000;'>                              tutorial9.py</span></i>
<i><span style='color:#008000;'>                 DAE Tools: pyDAE module, www.daetools.com</span></i>
<i><span style='color:#008000;'>                 Copyright (C) Dragan Nikolic, 2010</span></i>
<i><span style='color:#008000;'>***********************************************************************************</span></i>
<i><span style='color:#008000;'>DAE Tools is free software; you can redistribute it and/or modify it under the </span></i>
<i><span style='color:#008000;'>terms of the GNU General Public License as published by the Free Software </span></i>
<i><span style='color:#008000;'>Foundation; either version 3 of the License, or (at your option) any later version.</span></i>
<i><span style='color:#008000;'>The DAE Tools is distributed in the hope that it will be useful, but WITHOUT ANY </span></i>
<i><span style='color:#008000;'>WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A </span></i>
<i><span style='color:#008000;'>PARTICULAR PURPOSE. See the GNU General Public License for more details.</span></i>
<i><span style='color:#008000;'>You should have received a copy of the GNU General Public License along with the</span></i>
<i><span style='color:#008000;'>DAE Tools software; if not, write to the Free Software Foundation, Inc., </span></i>
<i><span style='color:#008000;'>59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.</span></i>
<i><span style='color:#008000;'>********************************************************************************&quot;&quot;&quot;</span></i>

<i><span style='color:#008000;'>&quot;&quot;&quot;</span></i>
<i><span style='color:#008000;'>In this example we use the same conduction problem as in the tutorial 1.</span></i>
<i><span style='color:#008000;'>Here we introduce:</span></i>
<i><span style='color:#008000;'> - Third party linear equations solvers</span></i>

<i><span style='color:#008000;'>Currently there are 3rd party linear equations solvers:</span></i>
<i><span style='color:#008000;'> - TrilinosAmesos: sequential sparse direct solver defined in pyTrilinosAmesos module (GNU Lesser GPL)</span></i>
<i><span style='color:#008000;'> - IntelPardiso: multi-threaded sparse direct solver defined in pyIntelPardiso module (proprietary)</span></i>
<i><span style='color:#008000;'> - AmdACML: multi-threaded dense lapack direct solver defined in pyAmdACML (proprietary)</span></i>
<i><span style='color:#008000;'> - IntelMKL: multi-threaded dense lapack direct solver defined in pyIntelMKL (proprietary)</span></i>
<i><span style='color:#008000;'> - Lapack: generic sequential dense lapack direct solver defined in pyLapack module </span></i>
<i><span style='color:#008000;'>           (The University of Tennessee free license)</span></i>
<i><span style='color:#008000;'> - Atlas: Automatically Tuned Linear Algebra Software implementation of the sequential dense Lapack </span></i>
<i><span style='color:#008000;'>          direct solver defined in pyAtlas module (BSD-style license)</span></i>
<i><span style='color:#008000;'>&quot;&quot;&quot;</span></i>

<span style='color:#0000ff;'>import</span> sys
<span style='color:#0000ff;'>from</span> daetools.pyDAE <span style='color:#0000ff;'>import</span> <span style='color:#0000ff;'>*</span>
<span style='color:#0000ff;'>from</span> time <span style='color:#0000ff;'>import</span> localtime, strftime

<i><span style='color:#008000;'># First import desired solver's module:</span></i>
<i><span style='color:#008000;'>#import daetools.pyTrilinosAmesos as pyTrilinosAmesos</span></i>
<i><span style='color:#008000;'>#import daetools.pyIntelPardiso   as pyIntelPardiso</span></i>
<i><span style='color:#008000;'>#import daetools.pyAmdACML        as pyAmdACML</span></i>
<i><span style='color:#008000;'>#import daetools.pyIntelMKL       as pyIntelMKL</span></i>
<i><span style='color:#008000;'>#import daetools.pyLapack         as pyLapack</span></i>
<i><span style='color:#008000;'>#import daetools.pyAtlas          as pyAtlas</span></i>

typeNone         <span style='color:#0000ff;'>=</span> daeVariableType(<span style='color:#bf0303;'>&quot;None&quot;</span>,         <span style='color:#bf0303;'>&quot;-&quot;</span>,      <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>1E10</span>,   <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>1e-5</span>)
typeTemperature  <span style='color:#0000ff;'>=</span> daeVariableType(<span style='color:#bf0303;'>&quot;Temperature&quot;</span>,  <span style='color:#bf0303;'>&quot;K&quot;</span>,    <span style='color:#c000c0;'>100</span>, <span style='color:#c000c0;'>1000</span>, <span style='color:#c000c0;'>300</span>, <span style='color:#c000c0;'>1e-5</span>)
typeConductivity <span style='color:#0000ff;'>=</span> daeVariableType(<span style='color:#bf0303;'>&quot;Conductivity&quot;</span>, <span style='color:#bf0303;'>&quot;W/mK&quot;</span>,   <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>1E10</span>, <span style='color:#c000c0;'>100</span>, <span style='color:#c000c0;'>1e-5</span>)
typeDensity      <span style='color:#0000ff;'>=</span> daeVariableType(<span style='color:#bf0303;'>&quot;Density&quot;</span>,      <span style='color:#bf0303;'>&quot;kg/m3&quot;</span>,  <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>1E10</span>, <span style='color:#c000c0;'>100</span>, <span style='color:#c000c0;'>1e-5</span>)
typeHeatCapacity <span style='color:#0000ff;'>=</span> daeVariableType(<span style='color:#bf0303;'>&quot;HeatCapacity&quot;</span>, <span style='color:#bf0303;'>&quot;J/KgK&quot;</span>,  <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>1E10</span>, <span style='color:#c000c0;'>100</span>, <span style='color:#c000c0;'>1e-5</span>)

<b>class</b> modTutorial(daeModel):
    <b>def</b> <b><span style='color:#000e52;'>__init__</span></b>(<span style='color:#0000ff;'>self</span>, Name, Parent <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>None</span>, Description <span style='color:#0000ff;'>=</span> <span style='color:#bf0303;'>&quot;&quot;</span>):
        daeModel.<b><span style='color:#000e52;'>__init__</span></b>(<span style='color:#0000ff;'>self</span>, Name, Parent, Description)

        <span style='color:#0000ff;'>self</span>.x  <span style='color:#0000ff;'>=</span> daeDomain(<span style='color:#bf0303;'>&quot;x&quot;</span>, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;X axis domain&quot;</span>)
        <span style='color:#0000ff;'>self</span>.y  <span style='color:#0000ff;'>=</span> daeDomain(<span style='color:#bf0303;'>&quot;y&quot;</span>, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Y axis domain&quot;</span>)

        <span style='color:#0000ff;'>self</span>.Qb <span style='color:#0000ff;'>=</span> daeParameter(<span style='color:#bf0303;'>&quot;Q_b&quot;</span>, eReal, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Heat flux at the bottom edge of the plate, W/m2&quot;</span>)
        <span style='color:#0000ff;'>self</span>.Qt <span style='color:#0000ff;'>=</span> daeParameter(<span style='color:#bf0303;'>&quot;Q_t&quot;</span>, eReal, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Heat flux at the top edge of the plate, W/m2&quot;</span>)

        <span style='color:#0000ff;'>self</span>.ro <span style='color:#0000ff;'>=</span> daeParameter(<span style='color:#bf0303;'>&quot;&amp;rho;&quot;</span>, eReal, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Density of the plate, kg/m3&quot;</span>)
        <span style='color:#0000ff;'>self</span>.cp <span style='color:#0000ff;'>=</span> daeParameter(<span style='color:#bf0303;'>&quot;c_p&quot;</span>, eReal, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Specific heat capacity of the plate, J/kgK&quot;</span>)
        <span style='color:#0000ff;'>self</span>.k  <span style='color:#0000ff;'>=</span> daeParameter(<span style='color:#bf0303;'>&quot;&amp;lambda;&quot;</span>,  eReal, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Thermal conductivity of the plate, W/mK&quot;</span>)
 
        <span style='color:#0000ff;'>self</span>.T <span style='color:#0000ff;'>=</span> daeVariable(<span style='color:#bf0303;'>&quot;T&quot;</span>, typeTemperature, <span style='color:#0000ff;'>self</span>, <span style='color:#bf0303;'>&quot;Temperature of the plate, K&quot;</span>)
        <span style='color:#0000ff;'>self</span>.T.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.x)
        <span style='color:#0000ff;'>self</span>.T.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.y)

    <b>def</b> DeclareEquations(<span style='color:#0000ff;'>self</span>):
        eq <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.CreateEquation(<span style='color:#bf0303;'>&quot;HeatBalance&quot;</span>, <span style='color:#bf0303;'>&quot;Heat balance equation. Valid on the open x and y domains&quot;</span>)
        x <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.x, eOpenOpen)
        y <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.y, eOpenOpen)
        eq.Residual <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.ro() <span style='color:#0000ff;'>*</span> <span style='color:#0000ff;'>self</span>.cp() <span style='color:#0000ff;'>*</span> <span style='color:#0000ff;'>self</span>.T.dt(x, y) <span style='color:#0000ff;'>-</span> <span style='color:#0000ff;'>self</span>.k() <span style='color:#0000ff;'>*</span> <span style='color:#0000ff;'>\</span>
                     (<span style='color:#0000ff;'>self</span>.T.d2(<span style='color:#0000ff;'>self</span>.x, x, y) <span style='color:#0000ff;'>+</span> <span style='color:#0000ff;'>self</span>.T.d2(<span style='color:#0000ff;'>self</span>.y, x, y))

        eq <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.CreateEquation(<span style='color:#bf0303;'>&quot;BC_bottom&quot;</span>, <span style='color:#bf0303;'>&quot;Boundary conditions for the bottom edge&quot;</span>)
        x <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.x, eClosedClosed)
        y <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.y, eLowerBound)
        eq.Residual <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>-</span> <span style='color:#0000ff;'>self</span>.k() <span style='color:#0000ff;'>*</span> <span style='color:#0000ff;'>self</span>.T.d(<span style='color:#0000ff;'>self</span>.y, x, y) <span style='color:#0000ff;'>-</span> <span style='color:#0000ff;'>self</span>.Qb()

        eq <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.CreateEquation(<span style='color:#bf0303;'>&quot;BC_top&quot;</span>, <span style='color:#bf0303;'>&quot;Boundary conditions for the top edge&quot;</span>)
        x <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.x, eClosedClosed)
        y <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.y, eUpperBound)
        eq.Residual <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>-</span> <span style='color:#0000ff;'>self</span>.k() <span style='color:#0000ff;'>*</span> <span style='color:#0000ff;'>self</span>.T.d(<span style='color:#0000ff;'>self</span>.y, x, y) <span style='color:#0000ff;'>-</span> <span style='color:#0000ff;'>self</span>.Qt()

        eq <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.CreateEquation(<span style='color:#bf0303;'>&quot;BC_left&quot;</span>, <span style='color:#bf0303;'>&quot;Boundary conditions at the left edge&quot;</span>)
        x <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.x, eLowerBound)
        y <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.y, eOpenOpen)
        eq.Residual <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.T.d(<span style='color:#0000ff;'>self</span>.x, x, y)

        eq <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.CreateEquation(<span style='color:#bf0303;'>&quot;BC_right&quot;</span>, <span style='color:#bf0303;'>&quot;Boundary conditions for the right edge&quot;</span>)
        x <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.x, eUpperBound)
        y <span style='color:#0000ff;'>=</span> eq.DistributeOnDomain(<span style='color:#0000ff;'>self</span>.y, eOpenOpen)
        eq.Residual <span style='color:#0000ff;'>=</span> <span style='color:#0000ff;'>self</span>.T.d(<span style='color:#0000ff;'>self</span>.x, x, y)

<b>class</b> simTutorial(daeDynamicSimulation):
    <b>def</b> <b><span style='color:#000e52;'>__init__</span></b>(<span style='color:#0000ff;'>self</span>):
        daeDynamicSimulation.<b><span style='color:#000e52;'>__init__</span></b>(<span style='color:#0000ff;'>self</span>)
        <span style='color:#0000ff;'>self</span>.m <span style='color:#0000ff;'>=</span> modTutorial(<span style='color:#bf0303;'>&quot;Tutorial_9&quot;</span>)
        <span style='color:#0000ff;'>self</span>.m.Description <span style='color:#0000ff;'>=</span> <span style='color:#bf0303;'>&quot;This tutorial explains how to create 3rd part linear solvers. &quot;</span>
          
    <b>def</b> SetUpParametersAndDomains(<span style='color:#0000ff;'>self</span>):
        n <span style='color:#0000ff;'>=</span> <span style='color:#c000c0;'>25</span>
        
        <span style='color:#0000ff;'>self</span>.m.x.CreateDistributed(eCFDM, <span style='color:#c000c0;'>2</span>, n, <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>0.1</span>)
        <span style='color:#0000ff;'>self</span>.m.y.CreateDistributed(eCFDM, <span style='color:#c000c0;'>2</span>, n, <span style='color:#c000c0;'>0</span>, <span style='color:#c000c0;'>0.1</span>)

        <span style='color:#0000ff;'>self</span>.m.k.SetValue(<span style='color:#c000c0;'>401</span>)
        <span style='color:#0000ff;'>self</span>.m.cp.SetValue(<span style='color:#c000c0;'>385</span>)
        <span style='color:#0000ff;'>self</span>.m.ro.SetValue(<span style='color:#c000c0;'>8960</span>)

        <span style='color:#0000ff;'>self</span>.m.Qb.SetValue(<span style='color:#c000c0;'>1e6</span>)
        <span style='color:#0000ff;'>self</span>.m.Qt.SetValue(<span style='color:#c000c0;'>0</span>)

    <b>def</b> SetUpVariables(<span style='color:#0000ff;'>self</span>):
        <b>for</b> x <span style='color:#0000ff;'>in</span> <b><span style='color:#000000;'>range</span></b>(<span style='color:#c000c0;'>1</span>, <span style='color:#0000ff;'>self</span>.m.x.NumberOfPoints <span style='color:#0000ff;'>-</span> <span style='color:#c000c0;'>1</span>):
            <b>for</b> y <span style='color:#0000ff;'>in</span> <b><span style='color:#000000;'>range</span></b>(<span style='color:#c000c0;'>1</span>, <span style='color:#0000ff;'>self</span>.m.y.NumberOfPoints <span style='color:#0000ff;'>-</span> <span style='color:#c000c0;'>1</span>):
                <span style='color:#0000ff;'>self</span>.m.T.SetInitialCondition(x, y, <span style='color:#c000c0;'>300</span>)

<b>if</b> <b><span style='color:#000000;'>__name__</span></b> <span style='color:#0000ff;'>==</span> <span style='color:#bf0303;'>&quot;__main__&quot;</span>:
    <i><span style='color:#008000;'># Create Log, Solver, DataReporter and Simulation object</span></i>
    log          <span style='color:#0000ff;'>=</span> daePythonStdOutLog()
    solver       <span style='color:#0000ff;'>=</span> daeIDASolver()
    datareporter <span style='color:#0000ff;'>=</span> daeTCPIPDataReporter()
    simulation   <span style='color:#0000ff;'>=</span> simTutorial()

    <i><span style='color:#008000;'># The default linear solver is Sundials dense sequential solver (LU decomposition).</span></i>
    <i><span style='color:#008000;'># It is possible to set the following 3rd party direct linear solvers:</span></i>
    <i><span style='color:#008000;'>#  1. Sparse solvers:</span></i>
    <i><span style='color:#008000;'>#      - IntelPardiso (multi-threaded - OMP)</span></i>
    <i><span style='color:#008000;'>#      - Trilinos Amesos (sequential): Klu, SuperLU, Lapack, Umfpack</span></i>
    <i><span style='color:#008000;'>#  3. Dense lapack wrappers:</span></i>
    <i><span style='color:#008000;'>#      - Amd ACML (OMP)</span></i>
    <i><span style='color:#008000;'>#      - Intel MKL (OMP)</span></i>
    <i><span style='color:#008000;'>#      - Generic Lapack (Sequential)</span></i>
    <i><span style='color:#008000;'>#      - Atlas lapack (Sequential)</span></i>
    <i><span style='color:#008000;'># If you are using Intel/AMD solvers you have to export their bin directories (see their docs how to do it).</span></i>
    <i><span style='color:#008000;'># If you are using OMP capable solvers you should set the number of threads to the number of cores. </span></i>
    <i><span style='color:#008000;'># For instance:</span></i>
    <i><span style='color:#008000;'>#    export OMP_NUM_THREADS=4</span></i>
    <i><span style='color:#008000;'># You can place the above command at the end of $HOME/.bashrc (or type it in shell, before simulation).</span></i>

    <i><span style='color:#008000;'># Import desired module and uncomment corresponding solver and set it by using SetLASolver function</span></i>
    <i><span style='color:#008000;'>#</span></i>
    <i><span style='color:#008000;'>#print &quot;Supported Trilinos Amesos 3rd party LA solvers:&quot;, str(pyTrilinosAmesos.daeTrilinosAmesosSupportedSolvers())</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyTrilinosAmesos.daeCreateTrilinosAmesosSolver(&quot;Amesos_Klu&quot;)</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyTrilinosAmesos.daeCreateTrilinosAmesosSolver(&quot;Amesos_Superlu&quot;)</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyTrilinosAmesos.daeCreateTrilinosAmesosSolver(&quot;Amesos_Lapack&quot;)</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyTrilinosAmesos.daeCreateTrilinosAmesosSolver(&quot;Amesos_Umfpack&quot;)</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyIntelPardiso.daeCreateIntelPardisoSolver()</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyAmdACML.daeCreateLapackSolver()</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyIntelMKL.daeCreateLapackSolver()</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyLapack.daeCreateLapackSolver()</span></i>
    <i><span style='color:#008000;'>#lasolver     = pyAtlas.daeCreateLapackSolver()</span></i>
    <i><span style='color:#008000;'>#solver.SetLASolver(lasolver)</span></i>

    <i><span style='color:#008000;'># Enable reporting of all variables</span></i>
    simulation.m.SetReportingOn(<span style='color:#0000ff;'>True</span>)

    <i><span style='color:#008000;'># Set the time horizon and the reporting interval</span></i>
    simulation.ReportingInterval <span style='color:#0000ff;'>=</span> <span style='color:#c000c0;'>10</span>
    simulation.TimeHorizon <span style='color:#0000ff;'>=</span> <span style='color:#c000c0;'>1000</span>

    <i><span style='color:#008000;'># Connect data reporter</span></i>
    simName <span style='color:#0000ff;'>=</span> simulation.m.Name <span style='color:#0000ff;'>+</span> strftime(<span style='color:#bf0303;'>&quot; [</span><span style='color:#0000ff;'>%d</span><span style='color:#bf0303;'>.%m.%Y %H:%M:%S]&quot;</span>, localtime())
    <b>if</b>(datareporter.Connect(<span style='color:#bf0303;'>&quot;&quot;</span>, simName) <span style='color:#0000ff;'>==</span> <span style='color:#0000ff;'>False</span>):
        sys.exit()

    <i><span style='color:#008000;'># Initialize the simulation</span></i>
    simulation.Initialize(solver, datareporter, log)

    <i><span style='color:#008000;'># Save the model report and the runtime model report </span></i>
    simulation.m.SaveModelReport(simulation.m.Name <span style='color:#0000ff;'>+</span> <span style='color:#bf0303;'>&quot;.xml&quot;</span>)
    simulation.m.SaveRuntimeModelReport(simulation.m.Name <span style='color:#0000ff;'>+</span> <span style='color:#bf0303;'>&quot;-rt.xml&quot;</span>)

    <i><span style='color:#008000;'># Solve at time=0 (initialization)</span></i>
    simulation.SolveInitial()

    <i><span style='color:#008000;'># Run</span></i>
    simulation.Run()
</pre>
</body>
</html>
